##########################
#Alternative Explanations#
##########################
#This to obtain:
#table_A18.tex

library("broom")
library("dplyr")
library("grid")
library(tidyr)
library(sandwich)
library(lfe)
library(broom)
library(dplyr)
library(tidyr)
library(grid)
library(stargazer)
library(lmtest)
library(multiwayvcov)
library("ggplot2")
library("glue")
library("spdep")
library("spatialreg")
library("stringr")
library("readxl")  
library("rgdal")

setwd("/Users/bgpopescu/")
#setwd("C:/Users/bogdanp/")

#############
#1921 Census#
#############

census_1921 <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                       layer="census_1921_points")

census_1921$dist1 <- as.numeric( with (census_1921,ifelse(census_1921$treat==1, 1, -1)))
census_1921$dist2 <-census_1921$dist1*(census_1921$krajna6_distance/1000)

census_1921$quad1<-ifelse((census_1921$lon > 18 & census_1921$lon< 20) & 
                            (census_1921$lat<47  & census_1921$lat>44), 1, 0)

census_1921$quad2<-ifelse((census_1921$lon > 16 & census_1921$lon< 18) & 
                            (census_1921$lat<47  & census_1921$lat>44), 1, 0)

census_1921$quad3<-ifelse((census_1921$lon > 14 & census_1921$lon< 16) & 
                            (census_1921$lat<47  & census_1921$lat>44), 1, 0)

census_1921$bfe1 <- ifelse(census_1921$krajna6_NEAR_FID == 1, 1,0)
census_1921$bfe2 <- ifelse(census_1921$krajna6_NEAR_FID == 2, 1,0)
census_1921$bfe3 <- ifelse(census_1921$krajna6_NEAR_FID == 3, 1,0)
census_1921$bfe4 <- ifelse(census_1921$krajna6_NEAR_FID == 4, 1,0)
census_1921$bfe5 <- ifelse(census_1921$krajna6_NEAR_FID == 5, 1,0)
census_1921$bfe6 <- ifelse(census_1921$krajna6_NEAR_FID == 6, 1,0)
census_1921$bfe7 <- ifelse(census_1921$krajna6_NEAR_FID == 7, 1,0)
census_1921$bfe8 <- ifelse(census_1921$krajna6_NEAR_FID == 8, 1,0)

census_1921$eth_frac_1921<-1-((census_1921$Serb_or_Croat/census_1921$Total)^2+
                           (census_1921$Slovenian/census_1921$Total)^2+
                           (census_1921$Czech/census_1921$Total)^2+
                           (census_1921$Rusyn/census_1921$Total)^2+
                           (census_1921$Pole/census_1921$Total)^2+
                           (census_1921$Russian/census_1921$Total)^2+
                           (census_1921$Hungarian/census_1921$Total)^2+
                           (census_1921$German/census_1921$Total)^2+
                           (census_1921$Albanian/census_1921$Total)^2+
                           (census_1921$Turk/census_1921$Total)^2+
                           (census_1921$Romanian/census_1921$Total)^2+
                           (census_1921$Italian/census_1921$Total)^2+
                           (census_1921$French/census_1921$Total)^2+
                           (census_1921$English/census_1921$Total)^2+
                           (census_1921$Other_and_unknown/census_1921$Total)^2)

                           
census_1921$POINT_X<-census_1921$lon
census_1921$POINT_Y<-census_1921$lat

census_1921$rel_frac_1921<-1-((census_1921$Christian_orthodox/census_1921$Total)^2+
(census_1921$Roman_catholic/census_1921$Total)^2+
(census_1921$Greek_catholic/census_1921$Total)^2+
(census_1921$Evangelic/census_1921$Total)^2+
(census_1921$Muslim/census_1921$Total)^2+
(census_1921$Jew/census_1921$Total)^2+
(census_1921$Other/census_1921$Total)^2+
(census_1921$no_confession_unknown/census_1921$Total)^2)
census_1921_data<-subset(census_1921, !is.na(treat))


#############
#1931 Census#
#############

census_1931 <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                       layer="census_1931_points")

census_1931$dist1 <- as.numeric( with (census_1931,ifelse(census_1931$treat==1, 1, -1)))
census_1931$dist2 <-census_1931$dist1*(census_1931$krajna6_distance/1000)

census_1931$quad1<-ifelse((census_1931$lon > 18 & census_1931$lon< 20) & 
                            (census_1931$lat<47  & census_1931$lat>44), 1, 0)

census_1931$quad2<-ifelse((census_1931$lon > 16 & census_1931$lon< 18) & 
                            (census_1931$lat<47  & census_1931$lat>44), 1, 0)

census_1931$quad3<-ifelse((census_1931$lon > 14 & census_1931$lon< 16) & 
                            (census_1931$lat<47  & census_1931$lat>44), 1, 0)

census_1931$bfe1 <- ifelse(census_1931$krajna6_NEAR_FID == 1, 1,0)
census_1931$bfe2 <- ifelse(census_1931$krajna6_NEAR_FID == 2, 1,0)
census_1931$bfe3 <- ifelse(census_1931$krajna6_NEAR_FID == 3, 1,0)
census_1931$bfe4 <- ifelse(census_1931$krajna6_NEAR_FID == 4, 1,0)
census_1931$bfe5 <- ifelse(census_1931$krajna6_NEAR_FID == 5, 1,0)
census_1931$bfe6 <- ifelse(census_1931$krajna6_NEAR_FID == 6, 1,0)
census_1931$bfe7 <- ifelse(census_1931$krajna6_NEAR_FID == 7, 1,0)
census_1931$bfe8 <- ifelse(census_1931$krajna6_NEAR_FID == 8, 1,0)


census_1931$eth_frac_1931<-1-((census_1931$Orthodox/census_1931$total)^2+
                              (census_1931$Roman_catholic/census_1931$total)^2+
                              (census_1931$Protestant/census_1931$total)^2+
                              (census_1931$Other_christian/census_1931$total)^2+
                              (census_1931$Muslim/census_1931$total)^2+
                              (census_1931$no_confession_or_unknown/census_1931$total)^2)

census_1931<-subset(census_1931, !is.na(treat))


###########
#2011 Data#
###########
#Reading Polygons

x<-st_layers(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb")

HRV_adm0 <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                    layer="HRV_adm0_wgs")

krajna <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                  layer="krajna_line_GCS")

HRV_adm3 <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                    layer="HRV_adm3")
HRV_adm3<-st_simplify(HRV_adm3, dTolerance = 0.005)
HRV_adm3<-subset(HRV_adm3, select = c(ID_2))

HRV_adm3_data = read_excel("./Dropbox/Legacies_Project/Analysis/data/merge.xlsx", sheet=1, col_names = TRUE, skip = 0)

final_sp<-left_join(HRV_adm3, HRV_adm3_data, by = c("ID_2"="ID_2"))
final_sp$bosnia_border_mun
data<-subset(final_sp, bosnia_border_mun==0)

data$treat <- as.numeric( with (data,ifelse((data$treat_KRAJ_kraj6==1), 1, 0)))
data$dist1 <- as.numeric( with (data,ifelse(data$treat==1, 1, -1)))
data$dist2 <-data$dist1*data$krajna6_distance

##########
#Analyses#
##########


rel_frac_1921<-lm(rel_frac_1921~treat + POINT_X + POINT_Y + zagreb_distance +
                    bfe1 + bfe2 + bfe3 + bfe4 + bfe5 + bfe6 + bfe7 + 
                    quad1 + quad2 + quad3, data = census_1921)


eth_frac_1931<-lm(eth_frac_1931~treat + POINT_X + POINT_Y + zagreb_distance +
                    bfe1 + bfe2 + bfe3 + bfe4 + bfe5 + bfe6 + bfe7 + 
                    quad1 + quad2 + quad3, data = census_1931)

eth_frac_1991<-lm(eth_frac_1991~treat + POINT_X + POINT_Y + zagreb_distance +
                    bfe1 + bfe2 + bfe3 + bfe4 + bfe5 + bfe6 + bfe7 + 
                    quad1 + quad2 + quad3, data = data)

eth_frac_2001<-lm(eth_frac_2001~treat + POINT_X + POINT_Y + zagreb_distance +
                    bfe1 + bfe2 + bfe3 + bfe4 + bfe5 + bfe6 + bfe7 + 
                    quad1 + quad2 + quad3, data = data)

eth_frac_2011<-lm(eth_frac_2001~treat + POINT_X + POINT_Y + zagreb_distance +
                    bfe1 + bfe2 + bfe3 + bfe4 + bfe5 + bfe6 + bfe7 + 
                    quad1 + quad2 + quad3, data = data)

transparency<-lm(transparency~treat + POINT_X + POINT_Y + zagreb_distance +
                    bfe1 + bfe2 + bfe3 + bfe4 + bfe5 + bfe6 + bfe7 + 
                    quad1 + quad2 + quad3, data = data)

######################
#Step1: RELEVANT VARS#
######################

relevant_vars<-c("rel_frac_1921",
                 "eth_frac_1931",
                 "eth_frac_1991",
                 "eth_frac_2001",
                 "eth_frac_2011",
                 "transparency")

#####################
#Step2: MEANS AND SD#
#####################


mean_rel_frac_1921<-round(mean(census_1921$rel_frac_1921, na.rm=T), digits = 3)
mean_eth_frac_1931<-round(mean(census_1931$eth_frac_1931, na.rm=T), digits = 3)
mean_eth_frac_1991<-round(mean(data$eth_frac_1991, na.rm=T), digits = 3)
mean_eth_frac_2001<-round(mean(data$eth_frac_2001, na.rm=T), digits = 3)
mean_eth_frac_2011<-round(mean(data$eth_frac_2011, na.rm=T), digits = 3)
mean_transparency<-round(mean(data$transparency, na.rm=T), digits = 3)
means<-c(mean_rel_frac_1921,
            mean_eth_frac_1931,
            mean_eth_frac_1991,
            mean_eth_frac_2001,
            mean_eth_frac_2011,
            mean_transparency)

sd_rel_frac_1921<-round(sd(census_1921$rel_frac_1921, na.rm=T), digits = 3)
sd_eth_frac_1931<-round(sd(census_1931$eth_frac_1931, na.rm=T), digits = 3)
sd_eth_frac_1991<-round(sd(data$eth_frac_1991, na.rm=T), digits = 3)
sd_eth_frac_2001<-round(sd(data$eth_frac_2001, na.rm=T), digits = 3)
sd_eth_frac_2011<-round(sd(data$eth_frac_2011, na.rm=T), digits = 3)
sd_transparency<-round(sd(data$transparency, na.rm=T), digits = 3)

sds<-c(sd_rel_frac_1921,
            sd_eth_frac_1931,
            sd_eth_frac_1991,
            sd_eth_frac_2001,
            sd_eth_frac_2011,
            sd_transparency)


####################
#Step3: LIST MODELS#
####################

list_models<-list(rel_frac_1921,
                  eth_frac_1931,
                  eth_frac_1991,
                  eth_frac_2001,
                  eth_frac_2011,
                  transparency)


######################
#Step5: Column labels#
######################


column_labels= c("\\shortstack{Religious\\\\ Fractionalization\\\\ 1921}",
                 "\\shortstack{Ethnic\\\\ Fractionalization\\\\ 1931}",
                 "\\shortstack{Ethnic\\\\ Fractionalization\\\\ 1991}",
                 "\\shortstack{Ethnic\\\\ Fractionalization\\\\ 2001}",
                 "\\shortstack{Ethnic\\\\ Fractionalization\\\\ 2011}",
                 "\\shortstack{Local\\\\ Government\\\\ Transparency, 2016}")


rel_frac=stargazer(list_models,
                    column.labels= column_labels,
                    keep = c("treat"),
                    se = list(coeftest(rel_frac_1921, cluster.vcov(rel_frac_1921, census_1921$Kotar_Municipality))[, 2],
                              coeftest(eth_frac_1931, cluster.vcov(eth_frac_1931, census_1931$srez_name))[, 2],
                              coeftest(eth_frac_1991, cluster.vcov(eth_frac_1991, data$NAME_1))[, 2],
                              coeftest(eth_frac_2001, cluster.vcov(eth_frac_1991, data$NAME_1))[, 2],
                              coeftest(eth_frac_2011, cluster.vcov(eth_frac_2011, data$NAME_1))[, 2],
                              coeftest(transparency, cluster.vcov(transparency, data$NAME_1))[, 2]),
                    covariate.labels=c("Military Territory"),
                    dep.var.caption = "Dependent variable",
                    dep.var.labels.include = F,
                    omit.stat=c("LL","ser","f", "rsq"), no.space=TRUE, float=FALSE,
                    add.lines=list(c("Mean", means),
                                   c("SD", sds),
                                   c("Boundary FE", rep("Yes", 7)),
                                   c("Region FE", rep("Yes", 7))))
note_latex <- "\\multicolumn{7}{l} {\\parbox[t]{25cm}{ \\ \\textit{Notes:} 
Coefficients and robust standard errors in parantheses. $^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01.
The units of analysis are settlements for the corresponding censuses. 
See codebook in the appendix for data sources and how the variables were calculated.}}"
rel_frac[grepl("Note",rel_frac)] <- note_latex
cat(rel_frac, file="./Dropbox/Legacies_Project/Paper/tables/table_A18.tex", sep="\n")
